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Abstract 

I study the spreading of infectious diseases on heterogeneous populations. I represent the pop- 
ulation structure by a contact-graph where vertices represent agents and edges represent disease 
transmission channels among them. The population heterogeneity is taken into account by the 
agent's subdivision in types and the mixing matrix among them. I introduce a type-network repre- 
sentation for the mixing matrix allowing an intuitive understanding of the mixing patterns and the 
analytical calculations. Using an iterative approach I obtain recursive equations for the probability 
distribution of the outbreak size as a function of time. I demonstrate that the expected outbreak 
size and its progression in time are determined by the largest eigenvalue of the reproductive num- 
ber matrix and the characteristic distance between agents on the contact-graph. Finally, I discuss 
the impact of intervention strategies to halt epidemic outbreaks. This work provides both a qual- 
itative understanding and tools to obtain quantitative predictions for the spreading dynamics on 
heterogeneous populations. 
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I. INTRODUCTION 



The globalization of human interactions have created a fertile ground for the fast and 
broad spread of infectious diseases, potentially leading to worldwide epidemics. We are thus 
force to understand the spreading of infectious diseases within this global scenario. Yet, 
the study of worldwide epidemics is challenging given the heterogeneity of the populations 
involved 

The first sign of heterogeneity is given by the variability of the reproductive number 
within or across populations [(J Q, Isj]. The reproductive number is defined as the number 
of secondary cases generated by a primary infected case within a population of susceptible 
individuals. In the case of sexually transmitted diseases the reproductive number is pro- 
portional to the rate of sexual partner acquisition and it exhibits wide fluctuations 
l[ 0, 11, Q|. In network based approaches the reproductive number is proportional 



to the node's degree 



13. 



14j and it exhibits wide fluctuations as well In the absence 



of biases among the connections between agents this heterogeneity is completely taken into 
account by the reproductive number distribution ^| . 

There are other properties beyond the reproductive number requiring the subdivision 
of a population in different classes or types. This includes but is not limited to age, geo- 
graphical location, social status and sexual behavior. In general these heterogeneities cannot 
be characterized by a single probability distribution. They require a multi-type approach 
with probability distributions characterizing each type and a mixing matrix describing the 
patterns of transmission among them. 

Multi-type models are difficult to deal with and are generally tackled using multi-agent 
simulations j^, U, 0, [3, 17, 1^- The advantage of multi-agent simulations is that we can 
consider several details and study their impact on the spreading dynamics. On the other 
hand, given the large number of variables and model parameters it is difficult to understand 
which are the key variables driving the system's dynamics. Therefore, analytical calculations 
are required to funnel the multi-agent simulations into specific regions of the parameters 
space. 

In this work I study the spreading of infectious diseases on multi-type networks. I take as 
starting point the static problem formation developed by Newman Q and the theory of 
age-dependent multi-type branching process j2yj - I develop these mathematical approaches 
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to accommodate some distinctive properties of real networks that have not previously con- 
sidered. In section |H] I introduce the basic framework. Focusing on the population structure 
I consider the contact-graph characterizing the detailed interactions among agents and, at a 
metapopulation level, the type-network characterizing the interactions among agent's types. 
Through some simple examples I illustrate the properties of the mixing matrix and its type- 
network representation. This section ends defining a branching process modeling a spanning 
tree from an index agent to all other agents in the contact-graph. In section ITTT1 1 characterize 
the local spreading dynamics from an agent to its contacts, taking the susceptible, infected, 
and removed (SIR) model study. Bringing together the underlying network struc- 

ture and the local transmission dynamics in section 11111 1 define a branching process that 
models the disease spreading dynamics. In section I1VI I extend the iterative approach for 



a single type 



21 



22, 



23J to accommodate the particularities of the multi-type case. Fo- 



cusing on the expected behavior, in section IS I obtain general equations determining the 
progression of the expected number of cumulative and new infections. Starting from these 
equations I analyze some limited cases. First, I derive the final expected outbreak size 
and, second, I analyze the time progression of the expected outbreak size for the case of a 
time homogeneous local transmission. In section IV CI I discuss the impact of the population 
heterogeneity on intervention strategies. I emphasize the role of the characteristic distance 
between agents to quantify the impact of intervention strategies on small-world populations. 
I also illustrate interventions targeting specific agent's types using a bipartite population as 
a case study. Finally, in section EU I provide an overview of the main results and discuss 
future directions. 



II. POPULATION STRUCTURE 



Consider a population of N agents that are susceptible to an infectious disease. By agent 
I mean any entity that could host and transmit the disease. Since we are interested on the 
transmission of infectious diseases among humans an agent is a human in the first place. 
For vector-borne diseases we could have in addition agents representing the intermediary 
host while for airborne diseases an agent could also represent a public place. The agents are 
assumed to be heterogeneous meaning that there are different agent classes or types according 
to their pattern of connectivity to other agents and/or to the speed at which they could 
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potentially transmit an infectious disease. For instance, human can be divided according to 
their age, social status and geographical location. Furthermore, in the case of vector- and 
air-born diseases there is an additional type given by the non-human intermediary. More 
precisely, let us assume that the agent population is divided in M types and there are N a 
agents of type a = 1, . . . , M, satisfying the normalization condition 

M 
a=l 

Note that within this work I use the indexes a,b, . . . for the agent's type. In the following 
I introduce two representations of the population structure at the agent and type levels, 
respectively. 



A. Contact-graph 

The contact- grap h takes precisely into account who could potentially transmit the disease 



to whom 



More precisely, 



Definition II. 1. The contact-graph is a labeled graph where vertices represent agents, 
edges represent the potential disease transmission channels among them, and the vertices 
are labeled according to the agent's type. 

The contact-graph represents the population mixing at the agent's level. Since there is a 
one-to-one relation between vertices and the corresponding agents I use these two terms 
interchangeable. 

All the information necessary to characterize a given graph is provided by its adjacency 
matrix. Yet, we should take into account the large size of real populations and their change 
in time. In general, the only way to achieve such a detailed description relies on agent- 
based simulations. My scope is to bypass this detailed description and focus on statistical 
properties that does not depend on the population structure details or their change in time. 
Yet, to achieve that I need to specify the time scale where these statistical properties are 
measured. 

Excluding the effect of patient isolation or any other intervention, the time scale that 
matters is the time interval from the infection of an agent to its death or recovery, i.e. 
the disease life time within an agent. At this point I intentionally exclude the effect of 



interventions, such as patient isolation, in order to achieve a more general approach. Their 
influence is taken into account when defining the disease spreading dynamics (see section lTTT|) . 
It is also worth mentioning that the disease life time is a random variable. Therefore, the 
statistical properties introduced below are the expectation after averaging over the disease 
life time distribution. 

The degree of a node is the total number of edges emanating from it regardless the type 
of the node at the other end. Let p£' be probability distribution that a type a node has 
degree k and denote by 



(k) a = Y.Pt ] k. (2) 



k=l 



its mean. Note that by allowing k to take values larger than one we are already taking into 



account the existence of concurrency 



3, 29, a3- 



To characterize the spreading process it is also relevant to determine the same distribution 
but for a vertex found and the end of an edge selected at random. This sampling introduces 
a bias towards nodes with higher degree resulting in the probability distribution 



(a) _ kp ( ^ 



with average excess degree 



(*)<~) . (4) 

k 

where the minus one subtracts the edge from where the node was reached. Associated with 
these two probability distributions we introduce the generating functions 

oo 

tfa(*) = l>JV, (5) 
fc=0 

oo 

K(x) = £<z£V- 1 . (6) 

k=l 

From the derivatives of U a (x) and V a (x) we obtain the moments of and respectively. 
For instance 



Ua(l) = (k) a , 
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(7) 



K(l) = (k)^ . (8) 

Since the agent population is finite there is a typical distance D between every two agents 
on the contact-graph. Social experiments such as the Kevin Bacon and Erdos numbers [3ll 

n 

or the Milgram experiment [32| reveal that social actors are separated by a small number of 
acquaintances ( "small- world" property j^). This observation is supported by theoretical 



34. 
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approaches demonstrating that D grows at most as log N in random graphs 
More recently it has been shown that for several real networks D actually decreases or 
remains constant as the network evolve and increases its size Thus, I explicitly take 
into account that D is finite. 

Example II. 2 (Poisson contact process). Let us assume that type a agents establish 
connections with other agents at a constant rate A a and that the disease life time is constant 
and equal to T. In this case we obtain a Poisson distribution for the agent's degree 

Furthermore, q£ ] = p< o) , (k) a = (k)t cess) = A a T, and U{x) = V(x) = e^ 1 ^ . 



B. Type- network 

At the metapopulation level the population structure of the is determined by the mixing 
patterns among the different agent's types. Given a type a agent and one of its edges let 
e a b be the probability that the agent at the other end is of type b ( mixing matrix). From 
the mixing matrix we can construct the type-network characterizing the metapopulation 
structure. 

Definition II. 3. Type-network: In the type-network a node represents a type, an arc is 
drawn from type a to b if e a b > 0, and the arc's weights are given by e a b- 

Note that since e aa may be nonzero the type- network may contain loops. Fig. Q shows some 
simple type-networks. The single-type case is represented by a node with a loop (Fig. QJi). 
A bipartite population is represented by two nodes with an incoming and an outgoing arc 
(Fig. [TJo). This example could model a heterosexual population with no other distinction 
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FIG. 1: Type-network representation of simple mixing matrices, a) Single-type population, b) 
Bipartite population, c) Fully mixed population with three types, c) Two cities (open circles) and 
the commuters among them (solid circle). The continuous/dashed lines represent intra/inter city 
connections. 

than gender or a metapopulation given by people and public places A fully mixed 

population is represented by a complete network (Fig. ^p). A less intuitive example is the 
type-network shown in Fig. [TJi, representing a population divided in two cities and the 
commuters between them. 



C. Annealed spanning tree 

Given a contact-graph, let us consider an epidemic outbreak starting from a single agent 
(index case). In the worth case scenario the disease propagates to all the agents that could be 
reached from the index case using the network connections. Thus, the outbreak is represented 
by a spanning or causal tree from the index case to all reachable agents. On this tree, the 
generation of an agent corresponds with the topological or hopping distance from the index 
case. This picture motivates the introduction of the following branching process: 

Definition II. 4. Multi-type Annealed Spanning Tree (AST) 

Consider a labeled contact-graph characterized by {N a , p£ } and the type-network {e a b}- 
The multi-type annealed spanning tree (AST) is the branching process satisfying the follow- 
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ing properties: 



1. The process start from an index case of type a G {1, . . . , M} at generation d = 0. The 
index case generates k sons with probability distribution p^ 1 . Each son is of type b 
with probability e ah . 

2. Each son at generation 1 < d < D generates k — 1 sons with probability distribution 

Each son is of type b with probability e a b. 

3. A son at generation d = D does not generate new sons. 



The term annealed means that we are not analyzing the true (quenched) spanning tree on 
the graph but a branching process with similar statistical properties. This approximation is 
particularly good if the contact-graph is continuously changing in time albeit the constancy 
of its statistical properties. A similar mathematical construction has been previously intro- 



duced by Newman 



19. 39]. The main difference here is the explicit consideration of the 



truncation distance D. Finally, it is worth noticing that all results derived below are exact 
for the multi-type AST but an approximation for the original population structure. 



III. SPREADING DYNAMICS 



To proceed further we should specify how the disease is transmitted from an agent to its 
neighbors in the contact-graph. Let r a b be the probability that an infected agent of type 
a infects a susceptible neighbor of type b. Within this work I assume that if e a b > then 
r a b > 0. Indeed, the absence of transmission between two types is taken into account by the 
corresponding matrix element of e. Upon infection we also need to specify when it takes 
place. Given a type a agent (primary case) and one of its neighbors of type b (secondary 
case), we define the generation time as the time elapse from the infection of the 

primary case to the infection of the secondary case provided it happens. I assume that the 
generation times are independent random variables with the distribution function 

G ab {r) = Prob (X^ < r) , (10) 
parameterized by the type of the primary and secondary cases. 



S 



Example III.l (SIR model). In the SIR model agents can be in the three exclusive 
states susceptible, infected and removed. A susceptible agent is one that have not become 
infected but it is susceptible to acquire the infection. An infected agent is one that have 
already acquired the disease and can potentially transmit the disease. A removed agent is 
one that has been previously infected but it is already excluded from the spreading process. 
Within this work the removal of an agent takes into account intervention strategies resulting 
in the isolation of infected individuals from the disease transmission chain. The death or 
"natural" recovery of infected agents was already taken into account during the definition 
of the contact-graph in subsection III Al 

Consider an agent i of type a and one of its neighbors j of type b. Let Y^'^ be the infection 
time of agent j by i in the absence of agent's removal and let G[ a,b \r) = Prob (y^^ < t\ 
be its distribution function. Furthermore, let zf^ be the removal time of agent % in the 
presence of agent's removal and let G^' b \r) = Prob ^Zj^ < r^j be its distribution function. 
The probability that agent j is infected by agent i by time t is given by 

b ab (t)= f i dG l {r)[l-G R {r)] . (11) 
Jo 

From this magnitude we obtain the probability that agent j gets infected by agent i no 
matter when 



r ab = lim b ah (t) . (12) 

t—>oo 

and the distribution of generation times 

G ab (r) = —b ab (r) . (13) 

Tab 

The SIR model could be further generalized taking immunization into account. In this 
case non-infected agents are divided into susceptible and immune. If s a is the probability 
that a type a agent is immune then the probability that agent j is infected by agent i by 
time t reads 

b ab (t) = (1 - s b ) f dG x {r) [1 - G R (r)} . (14) 
Jo 

Furthermore, the transmission probability r ab and the generation time distribution G ab (r) 
are obtained substituting this equation into (fT2j) and (JT3j) . respectively. 
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These examples illustrates how to calculate the transmission probability r a b and the 
generation time distribution G a b{r) from the standards models characterizing the spreading 
of infectious diseases. More important, by encapsulating the model details into r ab and 
G a f,(r) we can obtain general results that are independent of these details. Later on, we can 
analyze the particularities of each model. 



A. Multi-type age-dependent AST 

At this point the local spreading dynamics has been completely specified and we can 
super-impose it on the multi-type AST. 

Definition III. 2. Multi-type age-dependent AST 

The multi-type age-dependent AST is composed of two elements, a multi-type AST 111.41 
and a local spreading dynamics defined by {r a b, G a b(r)}. The global dynamics is then 
specified by the following rules 

1. The process starts with an infected agent of type a G {1, . . . , M} while all other agents 
are susceptible. 

2. An infected agent of type a infects each of its neighbors of type b with probability r a b 
and generation time distribution G a b(r). 



The age-dependent AST is a generalization of the Bellman-Harris |40J and Crum-Mode 



Jagers |4l|, |42j multi-type age-dependent branching processes. The key new element is the 
truncation at a maximum generation, allowing us to consider the small-world property of 
real networks. In spite of the similarities the mathematical framework I implement deviates 
substantially from these previous approaches. Indeed, I exploit this truncation making a 
backward iteration from the final generation D to the index case. 



IV. ITERATIVE APPROACH 



Consider a branch of the AST rooted on a type a agent, at generation d, that was infected 
at time zero. Let Pn d ' a ' b \t) be the probability distribution to find n infected type b agents 
at time t on that branch. In particular p^°' a ' b ^ (£) is the probability distribution of the total 
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number of infected type b agents at time t on the whole AST, given the index case was 
of type a. Based on the tree structure we can develop an iterative approach to compute 
p(d,a,b) ^ recurs ively. 

Lemma IV. 1. Consider a type a infected agent at generation d of the multi-type age- 
dependent AST. This agent has degree k with probability pj^ for d = and excess degree 
k — 1 with probability qjf* for < d < D. Let us index by a its neighbors on the next 
generation d + 1, where a G {1, . . . , k} for d = 0, a G {1, . . . , k — 1} for < d < D, and 
a G {0} for d = D. Then 



Pi°' a ' b) (t) = p { a) [5 ab 5 nl + (1 - S ab )S n0 ] 



k M 



+ Yl ■ • • Yl ^ELi n a +6 ab ,n II Yl e °* 

k=l ni=0 fife=0 a=l c=l 

x r ac [ dG ac {r)P^ b \t - r) + 8 na>0 [l - r ac G ac (t)} 



(15) 



P&**\t) = q{ a) [5 ab 8 nl + (l-5 ab )5, 



n0\ 



k-1 M 



+ Y < lk ) Y--- Yl 5 E^in Q+5ab ,n II Y e - 
=2 ni=0 rife_i=0 a=l c=l 

r ac f dG ac (r)Pl d a +1 > c > b \t - r) + 5 na , [l - r ac G ac (t)] 



(16) 



p(D,a,b) {t) = + (1 _ 5af))5n0 

Proof. Let n be the number of infected type 6 agents on a branch rooted at type a agent, 
and let n a be the infected type b agents on the branches rooted at each of its neighbors a. 
Then 

n = 5 a b + Y na ' ( 18 ) 

a 

where 5 ab takes into account if the root agent is or it is not of type b. The probability dis- 
tribution of n is given by the sum of all the possible combinations of the random variables 
n a that satisfy (fTHj) . Now, the root agent and its neighbors lie on a tree and therefore n a 
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are independent random variables. Furthermore, all agents at generation d+1 has the same 
statistical properties, i.e. n a are identically distributed random variables. Therefore, the 
probability of each combination is decomposed into the product of the probability distribu- 
tion of the number of infected agents of type b on the sub-branches rooted at each neighbor. 
Thus, taking into account that each neighbors is of type c with probability e ac we obtain 

P^ b \t) = p { a) [5 ab 5 nl + (1 - 5 ab )5 n0 ] 

oo oo oo k M 

k=l ni=0 rifc=0 a=l c=l 

P^ b \t) = q[ a) [S ab 6 nl + (1 - 5 ab )S nQ ] 

oo oo oo k— 1 M 

+ E E • ■ ■ E II E w . (M) 

fc=2 ni=0 nfe„!=0 a=l c=l 

where Qnl +1 '°' c ' (t) is the probability distribution of n a which we proceed to calculate. 

Let us focus on one neighbor and let us assume that it is of type c. With probability 
1 — r ac this agent is not infected at any time and with probability r ac [l — G ac {t)] it is not 
yet infected at time t given it will be infected at some later time, resulting in 

Q { o +1 ' a ' C ' b \t) = l-r ac G ac . (21) 

On the other hand, with probability r ac the neighbor is infected at some time r, with dis- 
tribution function G ac (r), and the spreading dynamics continue to subsequent generations. 
Once the neighbor is infected the number of infected agents of type b on that sub-branch is 
a random variable with probability distribution Pn d+1 ' c ' b \t — r). Therefore, for n > 

Q (A+l,a,c, b ) {f) = ^ t dGac{r)P {^ b) {t _ T y (22) 
JO 

Finally, substituting (|2*Tj) and (|2*2*j) into (JTHJ) and (|2~U|) we obtain equations (fTHjl and (JTBJ). 

The demonstration of (fTTj) is straightforward. For d = D the process stops and therefore 
there is only one infected agent, the root itself, which is or it is not of type b, resulting in 

□ 
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Associated with the probability distribution p^ d ' a ' b ^ (£) we introduce that generating func- 
tion 



F^ b \x,t) = ^P^ h \t) 



x ' 



(23) 



n=0 



Using the recursive relations for the probability distribution ([15j ) -([17 p we obtain the following 
recursive relations for the generating function 



/ m r t 
F^ h \x,t)=x 5 ^U a [Y^e ac l-r ac G ac (t)+r ac dG ac {r)F {l ^ h \x,t - 

\c=l L Jo 



(24) 



F^ h \x,t)=x^V a [y^e ac l-r ac G ac (t)+r ac / dG ac (r)F^ b \x,t - r) 

Vc=i L Jo 



(25) 



F {D ' a ' b \x,t) = x s « b . (26) 
These recursive equations are going to be useful in the following calculations. 



V. EXPECTED BEHAVIOR 



Given a infected agent of type a the expected number of secondary infections of type b 
it generates is given by 

Rab = (k) a e ab r ab (27) 

if it is the index case and by 

Rab = (k) { : xcess) e ab r ab (28) 

otherwise. The matrices R and R are extensions of the basic reproductive number to the 
multi-type case. In the following it becomes clear that R is more relevant and therefore I 
refer to it as the reproductive number matrix. 

Lemma V.l. Consider an ensemble of multi-type age-dependent AST \IIIJ% with index case 
of type a. Let N ab (t) be the mean total number of infected type b agents at time t and let 

13 



Iab(t)dt be the mean number of type b agents that are infected between time t and t + dt. 
Then 



D 



N ab {t)=Y,{H*r {d - 1) ) ab {t) , (29) 



d=l 
D 

(It 



a^El^^^.w. (3°) 



d=l 

where 



Habit) = RabGabif) , (31) 
Jab(t) = R ab G ab (t) , (32) 

and the multiplication symbolized by * involves a matrix multiplication and a convolution in 
time. For instance, 

(H * J) ab (t) = J2 drH ac {r)J cb {t - r) . (33) 

c=l Jo 

(r 2 ) ab {t) = {j*j) ab {t) . (34) 

Proof. Let 

N^ b \t) = gF(d ' a ' b)(1,t) (35) 
ox 

be the mean number of infected type b agents on the branch rooted at a type a agent 
at generation d. In particular, N ab (t) = N^ 0,a,b \t). Making use of the recursive relations 
(El- (El we obtain 



N^ a > b \t) = 5 ab + U a (l) r ac / dG ac (r)N^ 

c=l J( > 



(t - r) (36) 



N^(t)=5 ab + V a {l)Y^r ac / dG ac (r)N^ b \t-r) (37) 

c=l J 
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Iterating these recursive relations from d = D to d = we obtain (J29jl . Then differentiating 
with respect to time we finally obtain In this step we also make use of the relation 

between U(l) and V(l) and the average degrees (jZJ-flSJl- 

□ 

This Lemma provides explicit equations for the expected progression of an epidemic out- 
break. In some particular cases these equations may be further expressed in terms of elemen- 
tary functions allowing an straightforward interpretation. More generally these equations 
can be evaluated numerically in cases where further reduction is not possible. In addition, 
Theorem IV. H is a starting point for calculations addressing some limiting cases, which is the 
subject of the following subsections. 



A. Final outbreak size 



The final outbreak size is obtained taking the limit t — > oo in (|29|h resulting in 

D 



N ab (oc) = J2(RR d ~ 1 ) ab ■ (39) 



d=l 



When R can be diagonalized we can write R = PDF , where P is the matrix composed 
of the eigenvectors of R, D is the diagonal matrix constructed from the corresponding 
eigenvalues (p a , a = 1, . . . , M) and P -1 is the inverse of P. Thus is reduced to 

AUoo) = (RPNP- 1 ) , (40) 

V / ab 

where iV is a diagonal matrix with diagonal entries 

| fijzl f or Pa ^ 1 

N aa = I Pa - 1 ' T (41) 

[ D , for p a = 1 

The following two Theorems show that the only thing we need to estimate the order of 
magnitude of the expected outbreak size is the largest eigenvalue of the reproductive number 
matrix R. 
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Theorem V.2 (Complete type-network). Consider a complete type-network and let p 
be the largest eigenvalue of R Then 

o D - 1 

N ab (oo) = Uab^J—^ , (42) 

where u ab is indenpendent of D. 

Proof. The mixing matrix of a complete type-network is positive defined and, therefore, R 
(I27j) and R (I28|) are positive defined as well. From the Perron-Frobenius Theorem |43J] it 
follows that the largest eigenvalue of R is simple and all the entries of its corresponding left 
eigenvector v are different from zero and have the same sign. In particular we choose all the 
components of v to be positive such that 

M M „ 

(RR d ~ 1 ) ab = RaciR^U = ^ciR^U ■ (43) 

c=l c=l c 

Taking into account that Yl c v cRcb = pvb we obtain the inequalities 

(min) d-1 (r>r>d-l\ ^ (max) <2-i 

*o6 P < [RR Job < U ab p (44) 



where 



(min) • ^b 

i> ah = mmR ac — , 

c V n 



u (max) _ maXjRac !^ ) (45) 

C Vr 



From (J321) and (jHSJ) we obtain 



Finally, from this equation we obtain (|42|) with 



(45) 



(47) 



n . (min) , , (max) . / /0 \ 

< v> ab < u ab < u\ b < oo , (48) 
where the inequality u ab > follows from pHjl . 

□ 
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c" d 




FIG. 2: Strongly connected type-network with six types. The dashed lines indicate the possible 
paths from type a to b. Note that only the types c, d and c" are neighbors of type a and type b 
can be only reached from the last two. 

This result can be generalized to type-networks that may not be complete but are still 
strongly connected, i.e. there is a path from every type a to every type b. In this case 
some entries of R ac and (R d ~ l ) c b in (}4*3*)) may be zero. Intuitively this means that some 
types c may not be a neighbor of a and, if they are, there may not be a path from c to 
b (See Fig. |2J). More precisely, given a type a let Out(a) be its set of out-neighbors, i.e. 
Out(a) = {c\e ac > 0}, and given a type b let Ind{b) be the set of types from where b is 
reached after d hops on the type-network, i.e. In^ib) = {c\(e d ) c b > 0}. Furthermore, let 

S ( f' b) = Out{a) n In d ^{b) (49) 

denote the set of types that are out-neighbors of the index case type a and belong to at least 
one path of length d from a to b on the type-network. For instance, in the example in Fig. 
S[ a ' b) = 0, Si a ' b) = {c'}, Si a ' b) = {c"}, and S { d a ' b) ^ for all d > 3. 

Theorem V.3 (Strongly connected type-network). Consider a strongly connected type- 
network. Let p be the largest eigenvalue of R d a b the distance on the type-network from 
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type a to b, n = [D/d ab ] and D ab = nd ab . Then 



N ab (oo) = u ab ^ 



i-l 



(50) 



l<d<D\S ( . a ' b) yfH 



where u ab is independent of D. 



Proof. The conditions of the Perron- Frobenius theorem |43J are valid beyond positive defined 
matrices and holds for the mixing matrix representing a strongly connected network. Thus, 
the largest eigenvalue of R is simple and all the entries of its corresponding eigenvector v 
are different from zero and have the same sign. In particular we choose all the components 
of v to be positive. Based on this fact we can write (|43|) . There may be, however, some 
entries of e and thus of R (J27j) and that are zero. Indeed we can rewrite (f4"3"j) as 



(RR d 1 )ab 



E 



Rn 



-Vc(R d ~ 1 )cb ■ 



l<c<M|ceS^' 



(a,6) 



Thus (RR 



" x )a6 = whenever M = 0. Otherwise, we obtain the inequalities 



(51) 



where 



■ft 



(min) d—1 



ab 



p 4 - 1 < (RR 



d— In 



lab 



< U 



(max) d—1 



ab 



P 



(52) 



(min) . V b 

U ab = mm , s — R ac , 

ceOut(a) V c 



(53) 



From (|42j) and dHHJ) we obtain 



ft 



(max) 
ab 



max —R ac , 

ceOut(a) V c 



(54) 



1 + ^ E P"- 1 < N ab (oo) < 1 + u 

l<d<D\S^' b) ^% 

From this equation we obtain (|5()jl with 



(max) 
ab 



E 



l<d<D\S { d a ' b) ^tD 



(55) 



n . (min) ^ (max) 

< u\ b < u ab < u\ h < oo , 



where the inequality > follows from (f53l . 



(56) 



□ 
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FIG. 3: Expexted outbreak size as a function of the largest eigenvalue of the reproductive number 
matrix for different values of D. The region p < 1 is indicated by the dotted pattern. 

Figure El illustrates the predictions of Theorem 02] for complete type-networks. When 
p < 1 the expected outbreak size is of the order of the prefactor u a b which is not expected 
to be large. Different behaviors are observed, however, for p > 1 depending on D. For 
D 3> 1 there is a dramatic increase in the expected outbreak size. As soon as p > 1 a 
significant fraction of the agent population becomes affected. In contrast, when D is not so 
large it becomes clear that the expected outbreak size changes smoothly with increasing p, 
including the region around p — 1. This fact becomes relevant when analyzing the impact 
of intervention strategies (see section IV C|) . Finally, it is worth mentioning that a similar 
picture is obtained for the more general case of strongly-connected type-networks, albeit 
some corrections given by the missing terms the sum in (|50j). 

B. Spreading dynamics with constant transmission rate 



Now let us consider the particular case where the spreading dynamics is homogeneous, 
i.e. G a b(r) = G{r). In this case, from (j5Uj) we obtain the incidence 



D 

d=i 



(57) 
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In particular when R can be diagonalized we rewrite (|57j) as 



I ab (t) = (RPI(t)p- 

V / ab 

where I(t) is a time dependent diagonal matrix with diagonal entries 



(58) 



D 



(59) 



d=l 



Example V.4. Consider the case M = 2 with the reproductive number matrices 



R 



h 


k 2 




K x 


K 2 




, R = 




k 2 


h 




K 2 





(60) 



Since R is symmetric it can be diagonalized and P 1 = P T , where P T is the transpose of 
P. In this case R = PDP T with 



D 



Pi 

P 2 



P 



1 

71 



i i 
i -i 



(61) 



where 



p = Pl = K x + K 2 , p 2 = K 1 -K 2 (62) 
are the eigenvalues of R. Assuming an index case is of type a = 1 from (p)Hj) we finally obtain 

/ ll(f ) = *t±*"J u ( t ) + ^Ut) (63) 

h-Ai) = ^/„(t) - ^/ a {t) , (64) 

This example shows that in some cases we can exactly calculate the expected progression 
of an epidemic outbreak. More generally we obtain the following asymptotic behaviors. 

Theorem V.5. Consider a strongly connected type network and a homogeneous and expo- 
nential distribution of generation times G a b = 1 — e _Ar , where X is the transmission rate. 
Let p be the largest eigenvalue of R A2fy) and let 

(65) 
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9 > 1: If p > 1 and 1 < At < 9 then 



I ab (t) ~ e^ 1 ^ . (66) 

6> < 1: //At > 6» t/iera 



/ob(t) _ A(At) D -- 1 e" A ' 



iv a6 (oo) (D ab - 1; 

where D a b is the same as in Theorem \ V.3[ 



1 + 



(67) 



Proof. 9 ^> 1 : Following the same guidelines of the Theorem IV. 31 proof we arrive to the 
inequality 

u { r ] f ab (t) < U(t) < u^Ut) , (68) 

where 

/-co- E w 

l<d<_D|5^ a ' b) ^0 

The Laplace transform of f a b{t) is given by 

= r w*®*'** = - e (^a) d • ( ? °) 

P l<d<D\S ( f' b) ^H) 







When D — > oo this series converges only for > (p — 1)A. Therefore, / a b(t) ~ e^ -1 ^* when 
At — > oo. 

^ C 1: The demonstration of this case is straightforward. From Theorem IV. 31 it follows 
that (RR d ~ l )ab is of order p rf_1 for S'f'^ ^ 0. Therefore, forp^>D the sum in (JH7|) is 
dominated by the d = D a b term. Corrections are given by the ratio between the d = D a b 
and the preceeding term satisfying ^ 0, which is at most d = — 1. 

□ 

The case 9 ^> 1 provides the connection between this work and multi-type age-dependent 
branching processes with an infinite number of generations. Indeed, Mode have already 
demonstrated the exponential growth regime for the case D = oo (see j^j, Chapter 3). 
Theorem IV. 51 shows that on the other limit <C 1 the spreading dynamics is instead charac- 



terized by a gamma distribution, which is also the case for the single-type case 
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10 15 

t (days) 



FIG. 4: (a) Schematic representation of the evolution of a disease within an agent, starting from 
the moment the agent gets infected, passing through a latent state where the agent is not infectious 
and finally becoming infectious, (b) Gamma probability density function g(r) = G(t) for different 
values of a. (c) Number of secondary cases generated by a primary case for a SARS outbreak 
in Singapore, as reported in 44] (bars). The solid line is the best fit to the gamma probability 
density function times a pre-factor, resulting in a ~ 4.3. (d) Plot of the parameter 0(a) dividing 
the exponential and power law initial growth regimes as a function of the number of intermediary 
steps. 



Theorem IV.5I can be extended to consider other generation time distributions, such as a 
gamma distribution 



1 



Af 



r) = — - / dxx-'e-* , (71) 
r (a) Jo 

where a > 1. The gamma distribution can be interpreted as the existence of a — 1 inter- 
mediary steps before an agent becomes infectious (see Fig. EK,b). F° r en = 1 we recover 
the exponential distribution which corresponds with the absence of intermediary steps. The 
gamma distribution can be also obtained from the fit to some empirical distribution of 
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generation times (see Fig. Ht). 

In this case there are two important modifications to Theorem IV. 51 First, the parameter 
9 is now given by 

= [( nC -i)...( nJ - n) ]i 

P 

which increases approximately linearly with increasing a (see Fig. 0Ji). Second, in the regime 
6(a) 3> 1 although the fraction of infected agents is still given by a gamma distribution the 
exponent of the initial power law growth is given by aD, i.e. 

I ab (t) _\{\tr D ^e-^ 



N ab (<x>) T(aD ab - 1) 
Therefore, the existence of intermediary steps reduces the the small-world effect by a factor 
given by the number of intermediary steps a. For instance, by a factor of about four for 
SARS (Fig. Eb). 

C. Impact of intervention strategies 

The expected outbreak size is a monotonic increasing function of p (|50l) . which plays the 
role of the basic reproductive number in homogeneous populations [ll Lj]. Therefore, the 
aim of intervention strategies is to reduce the characteristic reproductive number p. On the 
other hand, intervention strategies implies an economical cost, including but not limited to 
the development of new vaccines and their deployment through vaccination campaigns. Our 
task is to design optimal intervention strategies that minimize the expected outbreak size 
with a feasible economical cost. 

To be more precise let us consider a scenario where the disease is transmitted at constant 
rate A from infected to susceptible agents, infected agents are isolated at a rate p and a 
fraction s of the population is immune to the disease. In this case the infection and removal 
times follows the exponential distribution functions G\{r) = 1 — e~ Ar and Ctr(t) = 1 — e~ MT , 
respectively. Thus, from (|12|) we obtain r ab = 1 — (3 where 

(3 = l--^(l-s), (74) 
\ + p 



23 



is the blocking fraction, i.e. the fraction of potential disease transmissions that are blocked 
either because of immunization or patient isolation. Since r ab = 1 — (3 is independent of the 
primary and secondary case types we can write the reproductive number matrices (J27j) and 
(125)) as R ab = (1 — (3)K and R — (1 — j3)K, respectively, where 

K ab = (k) a e ab , K ab = (k)^e ab . (75) 
In turn, the largest eigenvalue of R is given by 

P =(1-P)k, (76) 

where k is the largest eigenvalue of K. 

From the analysis made in section IV Al it follows that there are two different scenarios 
depending D ab . For simplicity let us focus on the complete type-network case where D ab = D. 
When D S> 1 the target of intervention strategies is p — 1, which is the consensus in the 



literature 



1Q- 



The blocking fraction to achieve this is obtained from 1)76)1. resulting in 

A = 1 - - ■ (77) 

n 

This result has been already reported, at least for the case of two types y|. When D is 
small, however, the expected outbreak size is a smooth function of p (see Fig. |3)). Therefore, 
(3 C does not represent a threshold value in small-world populations. 

So far we have considered homogenous intervention strategies. Now let us assume that 
the rate of patient isolation and the immunized fraction are now different for each agent's 
type and given by p a and s a , respectively. In this case the blocking fraction is given by 

/3 a6 = l-— ^— (l-s 6 ) , (78) 

A + Pa 

and r ab = 1 — /3 ab , which depends on the type of both the primary and secondary case. From 
the Perron- Frobenius Theorem it follows that p is a continuous increasing function of all the 

~ n 

entries of the corresponding matrix R [46]. Since R ab = (1 — f3 ab )K ab then p is a continuous 
decreasing function of 9 ab for all (a, b). The goal is to determine which strategy leads to the 
largest reduction of p. 

Example V.6. Consider the spread of HIV on an heterosexual population with no further 
distinction beyond gender. In this case the type-network is bipartite (see Fig. Wp)- Let k\ and 
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&2 be the average excess degree for the connections from women to men and biceversa. Let 
also assume that the rate of patient isolation is zero and that we could immunize a fraction 
s of the overall population, distributed between a fraction xs and (1 — x)s of immunized 
women and men, respectively. The question is to determine the value of x representing the 
best intervention strategy. In this case the reproductive number matrix is given by 



R 



and it has the largest eigenvalue 



P 



[l-(l-x)s)k 2 
(1 — xs)ki 



(79) 



— s + x(l — x)s 2 ]kik 2 . 



(80) 



It results that p is minimum for x = or x = 1, i.e. the best intervention strategy is to 
direct all the immunization resources to only one of the sub-populations. 

VI. DISCUSSION 



There is significant evidence that social networks are characterized by (i) wide connec- 
tivity fluctuations and (m) the small-world property 33]. The variability in the number of 
contacts (i) has a direct impact on the reproductive number. This fact has been taken into 
account since the seminal works of May and Anderson considering the variability in the rate 
of sexual partner acquisition^, 0, More recently it has gained attention for other infec- 
tious diseases as welL following the observation of super-spreading events in the 2002-2003 
SARS epidemics j?l all2|- Yet, the small- world property (ii) has been completely neglected. 



From my studies of the single type case 
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231 ] I have shown that intervention 



strategies are modulated by the average distance D between agents in the corresponding 
contact-graph. In this work I have demonstrated that this result is also valid for heteroge- 
neous populations. In this last case the characteristic reproductive number is given by the 
largest eigenvalue of the reproductive number matrix. The good news is that in spite of this 
modulation by D the target of intervention strategies is still the characteristic reproductive 
number. That is, the expected outbreak size still decreases with decreasing the characteris- 
tic reproductive number. The bad news is that to quantify the impact of the intervention 
strategies we need to estimate D. 
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There are different paths to estimate D. First, we can use a direct approach as the 
Milgram's experiments [32]. Second, we can measure other network proper t ies such as the 
degree distribution and then try to estimate D using network models 
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37 
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Finally, I have shown that the progression of the expected number of new infections is 



modulated by D (see 



22, 



23| and section UTT|) . More precissely, in small world populations 
the incidence is expected to grow as a power law and we can estimate D from the power law 
exponent. 

Further work is required to test the validity of the coarse grained description of the type- 
network approach. This can be done by running agent based simulations where we can have 
a strict control of the different statistical properties characterizing the population structure. 
These statistical properties can be then plug in into the type network approach to obtain 
qualitative and quantitative predictions that can be compared with the simulations results. 

In conclusion, this work opens new avenues to future research on the spreading of in- 
fectious diseases on heterogeneous populations. It allows for a qualitative understanding 
through the analysis of the type-network representation of the mixing matrix. More impor- 
tant, it leads to general results that can be tackled case by case using exact or approximate 
calculations and numerical computations. 

This work was supported by NSF Grants No. ITR 0426737 and No. ACT/SGER 0441089. 
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